package net.sourceforge.jaad.aac.sbr;

/**
 * This class is part of JAAD ( jaadec.sourceforge.net ) that is distributed
 * under the Public Domain license. Code changes provided by the JCodec project
 * are distributed under FreeBSD license.
 *
 * @author in-somnia
 */
class DCT {

	private static final int n = 32;

	// w_array_real[i] = cos(2*M_PI*i/32)
	private static final float[] w_array_real = {
		1.000000000000000f, 0.980785279337272f,
		0.923879528329380f, 0.831469603195765f,
		0.707106765732237f, 0.555570210304169f,
		0.382683402077046f, 0.195090284503576f,
		0.000000000000000f, -0.195090370246552f,
		-0.382683482845162f, -0.555570282993553f,
		-0.707106827549476f, -0.831469651765257f,
		-0.923879561784627f, -0.980785296392607f
	};

	// w_array_imag[i] = sin(-2*M_PI*i/32)
	private static final float[] w_array_imag = {
		0.000000000000000f, -0.195090327375064f,
		-0.382683442461104f, -0.555570246648862f,
		-0.707106796640858f, -0.831469627480512f,
		-0.923879545057005f, -0.980785287864940f,
		-1.000000000000000f, -0.980785270809601f,
		-0.923879511601754f, -0.831469578911016f,
		-0.707106734823616f, -0.555570173959476f,
		-0.382683361692986f, -0.195090241632088f
	};

	private static final float[] dct4_64_tab = {
		0.999924719333649f, 0.998118102550507f,
		0.993906974792480f, 0.987301409244537f,
		0.978317379951477f, 0.966976463794708f,
		0.953306019306183f, 0.937339007854462f,
		0.919113874435425f, 0.898674488067627f,
		0.876070082187653f, 0.851355195045471f,
		0.824589252471924f, 0.795836925506592f,
		0.765167236328125f, 0.732654273509979f,
		0.698376238346100f, 0.662415742874146f,
		0.624859452247620f, 0.585797846317291f,
		0.545324981212616f, 0.503538429737091f,
		0.460538715124130f, 0.416429549455643f,
		0.371317148208618f, 0.325310230255127f,
		0.278519600629807f, 0.231058135628700f,
		0.183039888739586f, 0.134580686688423f,
		0.085797272622585f, 0.036807164549828f,
		-1.012196302413940f, -1.059438824653626f,
		-1.104129195213318f, -1.146159529685974f,
		-1.185428738594055f, -1.221842169761658f,
		-1.255311965942383f, -1.285757660865784f,
		-1.313105940818787f, -1.337290763854981f,
		-1.358253836631775f, -1.375944852828980f,
		-1.390321016311646f, -1.401347875595093f,
		-1.408998727798462f, -1.413255214691162f,
		-1.414107084274292f, -1.411552190780640f,
		-1.405596733093262f, -1.396255016326904f,
		-1.383549690246582f, -1.367511272430420f,
		-1.348178386688232f, -1.325597524642944f,
		-1.299823284149170f, -1.270917654037476f,
		-1.238950133323669f, -1.203998088836670f,
		-1.166145324707031f, -1.125483393669128f,
		-1.082109928131104f, -1.036129593849182f,
		-0.987653195858002f, -0.936797380447388f,
		-0.883684754371643f, -0.828443288803101f,
		-0.771206021308899f, -0.712110757827759f,
		-0.651300072669983f, -0.588920354843140f,
		-0.525121808052063f, -0.460058242082596f,
		-0.393886327743530f, -0.326765477657318f,
		-0.258857429027557f, -0.190325915813446f,
		-0.121335685253143f, -0.052053272724152f,
		0.017354607582092f, 0.086720645427704f,
		0.155877828598022f, 0.224659323692322f,
		0.292899727821350f, 0.360434412956238f,
		0.427100926637650f, 0.492738455533981f,
		0.557188928127289f, 0.620297133922577f,
		0.681910991668701f, 0.741881847381592f,
		0.800065577030182f, 0.856321990489960f,
		0.910515367984772f, 0.962515234947205f,
		1.000000000000000f, 0.998795449733734f,
		0.995184719562531f, 0.989176511764526f,
		0.980785250663757f, 0.970031261444092f,
		0.956940352916718f, 0.941544055938721f,
		0.923879504203796f, 0.903989315032959f,
		0.881921231746674f, 0.857728600502014f,
		0.831469595432281f, 0.803207516670227f,
		0.773010432720184f, 0.740951120853424f,
		0.707106769084930f, 0.671558916568756f,
		0.634393274784088f, 0.595699310302734f,
		0.555570185184479f, 0.514102697372437f,
		0.471396654844284f, 0.427555114030838f,
		0.382683426141739f, 0.336889833211899f,
		0.290284633636475f, 0.242980122566223f,
		0.195090234279633f, 0.146730497479439f,
		0.098017133772373f, 0.049067649990320f,
		-1.000000000000000f, -1.047863125801086f,
		-1.093201875686646f, -1.135906934738159f,
		-1.175875544548035f, -1.213011503219605f,
		-1.247225046157837f, -1.278433918952942f,
		-1.306562900543213f, -1.331544399261475f,
		-1.353317975997925f, -1.371831417083740f,
		-1.387039899826050f, -1.398906826972961f,
		-1.407403707504273f, -1.412510156631470f,
		0f, -1.412510156631470f,
		-1.407403707504273f, -1.398906826972961f,
		-1.387039899826050f, -1.371831417083740f,
		-1.353317975997925f, -1.331544399261475f,
		-1.306562900543213f, -1.278433918952942f,
		-1.247225046157837f, -1.213011384010315f,
		-1.175875544548035f, -1.135907053947449f,
		-1.093201875686646f, -1.047863125801086f,
		-1.000000000000000f, -0.949727773666382f,
		-0.897167563438416f, -0.842446029186249f,
		-0.785694956779480f, -0.727051079273224f,
		-0.666655659675598f, -0.604654192924500f,
		-0.541196048259735f, -0.476434230804443f,
		-0.410524487495422f, -0.343625843524933f,
		-0.275899350643158f, -0.207508206367493f,
		-0.138617098331451f, -0.069392144680023f,
		0f, 0.069392263889313f,
		0.138617157936096f, 0.207508206367493f,
		0.275899469852448f, 0.343625962734222f,
		0.410524636507034f, 0.476434201002121f,
		0.541196107864380f, 0.604654192924500f,
		0.666655719280243f, 0.727051138877869f,
		0.785695075988770f, 0.842446029186249f,
		0.897167563438416f, 0.949727773666382f
	};

	private static final int[] bit_rev_tab = {0, 16, 8, 24, 4, 20, 12, 28, 2, 18, 10, 26, 6, 22, 14, 30, 1, 17, 9, 25, 5, 21, 13, 29, 3, 19, 11, 27, 7, 23, 15, 31};

	// FFT decimation in frequency
	// 4*16*2+16=128+16=144 multiplications
	// 6*16*2+10*8+4*16*2=192+80+128=400 additions
	private static void fft_dif(float[] Real, float[] Imag) {
		float w_real, w_imag; // For faster access
		float point1_real, point1_imag, point2_real, point2_imag; // For faster access
		int j, i, i2, w_index; // Counters

		// First 2 stages of 32 point FFT decimation in frequency
		// 4*16*2=64*2=128 multiplications
		// 6*16*2=96*2=192 additions
		// Stage 1 of 32 point FFT decimation in frequency
		for(i = 0; i<16; i++) {
			point1_real = Real[i];
			point1_imag = Imag[i];
			i2 = i+16;
			point2_real = Real[i2];
			point2_imag = Imag[i2];

			w_real = w_array_real[i];
			w_imag = w_array_imag[i];

			// temp1 = x[i] - x[i2]
			point1_real -= point2_real;
			point1_imag -= point2_imag;

			// x[i1] = x[i] + x[i2]
			Real[i] += point2_real;
			Imag[i] += point2_imag;

			// x[i2] = (x[i] - x[i2]) * w
			Real[i2] = ((point1_real*w_real)-(point1_imag*w_imag));
			Imag[i2] = ((point1_real*w_imag)+(point1_imag*w_real));
		}
		// Stage 2 of 32 point FFT decimation in frequency
		for(j = 0, w_index = 0; j<8; j++, w_index += 2) {
			w_real = w_array_real[w_index];
			w_imag = w_array_imag[w_index];

			i = j;
			point1_real = Real[i];
			point1_imag = Imag[i];
			i2 = i+8;
			point2_real = Real[i2];
			point2_imag = Imag[i2];

			// temp1 = x[i] - x[i2]
			point1_real -= point2_real;
			point1_imag -= point2_imag;

			// x[i1] = x[i] + x[i2]
			Real[i] += point2_real;
			Imag[i] += point2_imag;

			// x[i2] = (x[i] - x[i2]) * w
			Real[i2] = ((point1_real*w_real)-(point1_imag*w_imag));
			Imag[i2] = ((point1_real*w_imag)+(point1_imag*w_real));

			i = j+16;
			point1_real = Real[i];
			point1_imag = Imag[i];
			i2 = i+8;
			point2_real = Real[i2];
			point2_imag = Imag[i2];

			// temp1 = x[i] - x[i2]
			point1_real -= point2_real;
			point1_imag -= point2_imag;

			// x[i1] = x[i] + x[i2]
			Real[i] += point2_real;
			Imag[i] += point2_imag;

			// x[i2] = (x[i] - x[i2]) * w
			Real[i2] = ((point1_real*w_real)-(point1_imag*w_imag));
			Imag[i2] = ((point1_real*w_imag)+(point1_imag*w_real));
		}

		// Stage 3 of 32 point FFT decimation in frequency
		// 2*4*2=16 multiplications
		// 4*4*2+6*4*2=10*8=80 additions
		for(i = 0; i<n; i += 8) {
			i2 = i+4;
			point1_real = Real[i];
			point1_imag = Imag[i];

			point2_real = Real[i2];
			point2_imag = Imag[i2];

			// out[i1] = point1 + point2
			Real[i] += point2_real;
			Imag[i] += point2_imag;

			// out[i2] = point1 - point2
			Real[i2] = point1_real-point2_real;
			Imag[i2] = point1_imag-point2_imag;
		}
		w_real = w_array_real[4]; // = sqrt(2)/2
		// w_imag = -w_real; // = w_array_imag[4]; // = -sqrt(2)/2
		for(i = 1; i<n; i += 8) {
			i2 = i+4;
			point1_real = Real[i];
			point1_imag = Imag[i];

			point2_real = Real[i2];
			point2_imag = Imag[i2];

			// temp1 = x[i] - x[i2]
			point1_real -= point2_real;
			point1_imag -= point2_imag;

			// x[i1] = x[i] + x[i2]
			Real[i] += point2_real;
			Imag[i] += point2_imag;

			// x[i2] = (x[i] - x[i2]) * w
			Real[i2] = (point1_real+point1_imag)*w_real;
			Imag[i2] = (point1_imag-point1_real)*w_real;
		}
		for(i = 2; i<n; i += 8) {
			i2 = i+4;
			point1_real = Real[i];
			point1_imag = Imag[i];

			point2_real = Real[i2];
			point2_imag = Imag[i2];

			// x[i] = x[i] + x[i2]
			Real[i] += point2_real;
			Imag[i] += point2_imag;

			// x[i2] = (x[i] - x[i2]) * (-i)
			Real[i2] = point1_imag-point2_imag;
			Imag[i2] = point2_real-point1_real;
		}
		w_real = w_array_real[12]; // = -sqrt(2)/2
		// w_imag = w_real; // = w_array_imag[12]; // = -sqrt(2)/2
		for(i = 3; i<n; i += 8) {
			i2 = i+4;
			point1_real = Real[i];
			point1_imag = Imag[i];

			point2_real = Real[i2];
			point2_imag = Imag[i2];

			// temp1 = x[i] - x[i2]
			point1_real -= point2_real;
			point1_imag -= point2_imag;

			// x[i1] = x[i] + x[i2]
			Real[i] += point2_real;
			Imag[i] += point2_imag;

			// x[i2] = (x[i] - x[i2]) * w
			Real[i2] = (point1_real-point1_imag)*w_real;
			Imag[i2] = (point1_real+point1_imag)*w_real;
		}

		// Stage 4 of 32 point FFT decimation in frequency (no multiplications)
		// 16*4=64 additions
		for(i = 0; i<n; i += 4) {
			i2 = i+2;
			point1_real = Real[i];
			point1_imag = Imag[i];

			point2_real = Real[i2];
			point2_imag = Imag[i2];

			// x[i1] = x[i] + x[i2]
			Real[i] += point2_real;
			Imag[i] += point2_imag;

			// x[i2] = x[i] - x[i2]
			Real[i2] = point1_real-point2_real;
			Imag[i2] = point1_imag-point2_imag;
		}
		for(i = 1; i<n; i += 4) {
			i2 = i+2;
			point1_real = Real[i];
			point1_imag = Imag[i];

			point2_real = Real[i2];
			point2_imag = Imag[i2];

			// x[i] = x[i] + x[i2]
			Real[i] += point2_real;
			Imag[i] += point2_imag;

			// x[i2] = (x[i] - x[i2]) * (-i)
			Real[i2] = point1_imag-point2_imag;
			Imag[i2] = point2_real-point1_real;
		}

		// Stage 5 of 32 point FFT decimation in frequency (no multiplications)
		// 16*4=64 additions
		for(i = 0; i<n; i += 2) {
			i2 = i+1;
			point1_real = Real[i];
			point1_imag = Imag[i];

			point2_real = Real[i2];
			point2_imag = Imag[i2];

			// out[i1] = point1 + point2
			Real[i] += point2_real;
			Imag[i] += point2_imag;

			// out[i2] = point1 - point2
			Real[i2] = point1_real-point2_real;
			Imag[i2] = point1_imag-point2_imag;
		}

		//FFTReorder(Real, Imag);
	}

	/* size 64 only! */
	public static void dct4_kernel(float[] in_real, float[] in_imag, float[] out_real, float[] out_imag) {
		// Tables with bit reverse values for 5 bits, bit reverse of i at i-th position
		int i, i_rev;

		/* Step 2: modulate */
		// 3*32=96 multiplications
		// 3*32=96 additions
		for(i = 0; i<32; i++) {
			float x_re, x_im, tmp;
			x_re = in_real[i];
			x_im = in_imag[i];
			tmp = (x_re+x_im)*dct4_64_tab[i];
			in_real[i] = (x_im*dct4_64_tab[i+64])+tmp;
			in_imag[i] = (x_re*dct4_64_tab[i+32])+tmp;
		}

		/* Step 3: FFT, but with output in bit reverse order */
		fft_dif(in_real, in_imag);

		/* Step 4: modulate + bitreverse reordering */
		// 3*31+2=95 multiplications
		// 3*31+2=95 additions
		for(i = 0; i<16; i++) {
			float x_re, x_im, tmp;
			i_rev = bit_rev_tab[i];
			x_re = in_real[i_rev];
			x_im = in_imag[i_rev];

			tmp = (x_re+x_im)*dct4_64_tab[i+3*32];
			out_real[i] = (x_im*dct4_64_tab[i+5*32])+tmp;
			out_imag[i] = (x_re*dct4_64_tab[i+4*32])+tmp;
		}
		// i = 16, i_rev = 1 = rev(16);
		out_imag[16] = (in_imag[1]-in_real[1])*dct4_64_tab[16+3*32];
		out_real[16] = (in_real[1]+in_imag[1])*dct4_64_tab[16+3*32];
		for(i = 17; i<32; i++) {
			float x_re, x_im, tmp;
			i_rev = bit_rev_tab[i];
			x_re = in_real[i_rev];
			x_im = in_imag[i_rev];
			tmp = (x_re+x_im)*dct4_64_tab[i+3*32];
			out_real[i] = (x_im*dct4_64_tab[i+5*32])+tmp;
			out_imag[i] = (x_re*dct4_64_tab[i+4*32])+tmp;
		}

	}
}
